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J| BACKGROUND OF THE INVENTION 

§* Field of Invention 

The present invention relates generally to analyzing 

integrated circuits and, more particularly, to techniques for 
30 performing simulations of radio frequency (RF) integrated 

circuits . 

Discussion of Background 

The exploding demand for high performance wireless products 
35 has increased the need for efficient and accurate simulation 
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techniques for RF integrated circuits. RF circuit simulation is 
difficult because RF circuits typically contain signals with 
multiple-timescale properties, as usually the data and carrier 
signals in a system are separated in frequency by several orders 
of magnitude. Special-purpose RF simulators exploit the 

sparsity of the spectrum in order to make the computations 
tractable. See K. S. Kundert, "Introduction to RF Simulation 
and Its Applications," IEEE J. Sol. State Circuits, September 
1999. For example, instead of performing a long transient 
analysis of a circuit driven by a periodic source, we may seek 
to find the steady-state directly. Unfortunately, periodic- 
steady-state computation may be problematic in that it can be 
inefficient and inaccurate using conventional methods. 

Two numerical methods commonly used in steady-state 
computation are the shooting-Newton method, based on low-order 
finite difference discretizations such as the second-order Gear 
method (see K. S. Kundert, "Introduction to RF Simulation and 
Its Applications," IEEE J. Sol. State Circuits, September 1999), 
and the harmonic balance method, based on high-order spectral 
discretizations (see "Nonlinear Circuit Analysis Using The 
Method of Harmonic Balance-A Review of The Art," Part I— 
Introductory Concepts, Int. J. Microwave and Millimeter Wave 
Computer Aided Engineering, vol. 1, no. 1, 1991). 



An advantage of the low-order polynomial-based methods is 
that these methods can select time-points based on localized 
error estimates and as a result can easily handle sharp 
transitions in circuit waveforms. The harmonic balance method, 
on the other hand, has the advantage of attaining spectral 
accuracy for smooth waveforms. Recent developments in matrix- 
free Krylov-subspace algorithms have made these methods even 
more popular as they can be used to analyze circuits with 
thousands of devices. For more information on matrix-free 
Kylov-subspace algorithms, see R. Telichevesky, K. Kundert and 
J. K. White, "Efficient AC and Noise Analysis of Two-Tone RF 
Circuits, " Proc. 33rd Design Automation Conference, June, 1996; 
D. Long and R. Melville and K. Ashby and B. Horton, "Full Chip 
Harmonic Balance, " Proc. Custom Integrated Circuits Conference, 
May, 1997. 

High precision computations are often necessary in RF 
circuit simulation. For example, accurately computing the noise 
figure of a highly nonlinear RF circuit often requires accurate 
determination of the periodic operating point up to many 
multiples of the fundamental frequency, as noise may be 
translated from very high frequencies by the mixing action of 
the time-varying circuit elements to appear at the output. 

High accuracy is easy to achieve for smooth waveforms using 



spectral approximation in the harmonic balance method. However, 
non-linear circuits often produce waveforms that have sharp- 
transition regions. The waveforms may only be C° in the case of 
a pulse, which means that the waveforms are continuous but the 
derivatives of the waveforms are not continuous at some points. 
In this case, the spectral accuracy of the harmonic balance 
method is lost and acceptable accuracy is obtained only at the 
cost of including very many harmonics and/or timepoints in the 
simulation . 

To help handle sharp transitions, Nastov and White 
introduced the time-mapped harmonic balance (TMHB) method. See 
0. J. Nastov and J. K. White, "Time-Mapped Harmonic Balance," 
Proc. 36th Design Automation Conference, New Orleans, LA, June, 
1999. In the TMHB method, a time-map function is used to map 
the solution in the original time space to a space where the 
solution is more smooth. This method is an improvement over 
traditional harmonic balance, but is still restricted, since, 
theoretically speaking, a C° point in the original time space is 
still a C° point in the mapped time space, and so one still 
cannot achieve spectral, accuracy. Basically, the TMHB method 
can only put more points around the C° points to localize the 
oscillations inevitably associated with using a spectral method 
to approximate C° functions. And, as with all methods that 



utilize global basis functions, for reasons of numerical 
stability there are also strong restrictions on how rapidly the 
timestep can be changed locally. 

The low-order multi-step methods, such as the Gear methods, 
typically used in the shooting method offer excellent 
flexibility in locally adapting the timesteps. For low to 
moderate accuracy, these methods are preferred for problems that 
are highly nonlinear and/or contain sharp transitions. However, 
low-order multi-step methods are not efficient at high precision 
because very fine discretizations are needed, resulting in loss 
of speed and increased memory requirements. In addition, the 
multistep discretizations used in the shooting method must be 
causal, and so have stability problems associated with using 
one-side approximations, namely, performing differentiation 
using backward differences. It is generally accepted that 
multistep methods of order greater then three or four are not 
sufficiently stable, and in practice even the third order 
methods have relatively stringent stability restrictions on how 
rapidly timesteps may be varied. As a result, multi-step 
methods of order higher than 2 are not widely used in circuit 
simulation. Often, because of the backward-looking nature of 
the approximation, it is necessary to decrease the order after 
sharp-transition points or C° points, often to first order, 



resulting in very fine grids in those regions. 

In sum, most RF circuit analysis tools use either shooting- 
Newton or harmonic balance methods. Unfortunately , neither 
method can efficiently achieve high accuracy on strongly 
nonlinear circuits possessing waveforms with rapid transitions. 

SUMMARY OF THE INVENTION 

It has been recognized that what is needed is a method 
for perform efficient, highly accurate radio frequency circuit 
simulation. Broadly speaking, the present invention fills 
these needs by providing a method and device for multi- 
interval Chebyshev collocation. It should be appreciated that 
the present invention can be implemented in numerous ways, 
including as a process, an apparatus, a system, a device or a 
method. Several inventive embodiments of the present 

invention are described below. 

The present invention provides a multi-interval-Chebyshev 
(MIC) method for solving boundary-value differential-algebraic 
equations arising in radio frequency (RF) circuit simulation. 
The varying-order, varying-timestep MIC scheme effectively 
adapts to different types of waveforms that appear in circuit 
simulation. The method has spectral accuracy for smooth wave- 
forms. In some embodiments, the method is efficient in 
resolving sharp transitions and nonlinear effects with high 



accuracy. Note that purely low-order schemes have difficulty 
achieving accuracy at reasonable cost, and high-order methods 
such as harmonic balance lose accuracy due to oscillations in 
sharp-transition regions. Preconditioning techniques are 
presented herein that make the cost of the present scheme 
comparable to backward differential formula-based (BDF-based) 
shooting-Newton methods. 

In one embodiment, a method of simulating a circuit is 
provided. The method comprises the following: defining a 
differential-algebraic equation of the circuit; defining a 
simulation time interval corresponding to the differential- 
algebraic equation; dividing the simulation time interval into 
time intervals, wherein the time intervals have corresponding 
polynomials for each time interval, wherein each polynomial is 
a portion of an approximation to a desired solution of the 
differential-algebraic equation; and applying a collocation 
method to discretize the differential-algebraic equation. 

In another embodiment, a method of solving a set of 
differential-algebraic equation arising in a circuit 
simulation is provided. The method comprises the following: 
applying a collocation method to each differential-algebraic 
equation to discretize the set of differential-algebraic 
equations; and forming a solution to the set of differential- 



algebraic equations based on the discretized differential- 
algebraic equation. 

Advantageously , the present invention provides a method 
for efficient, high-accuracy radio frequency (RF) circuit 
simulation. Because of inevitable discontinuities in the 
solution waveforms, residuals, or device equations, the 
discretization is performed by breaking the simulation 
interval, of length T, into m timestep intervals of length h sr 
s = 1 . . . m. Within each interval, Chebyshev polynomials 
are used to represent the solution waveform. Chebyshev 
polynomials have excellent approximation properties, in 
particular showing spectral convergence as the order is 
increased. The method of the present invention has several 
desirable properties that lead to an excellent tradeoff 
between efficiency and accuracy. Further, the method can be 
easily combined with matrix-implicit Krylov-subspace 
techniques to allow analysis of large circuits. In a sense, 
the best characteristics of the shooting-Newton and harmonic- 
balance methods are exploited, without the drawbacks of 
either. 

BRIEF DESCRIPTION OF THE DRAWINGS 

The present invention will be readily understood by the 
following detailed description in conjunction with the 
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accompanying drawings. To facilitate this description, like 
reference numerals designate like structural elements. 

FIG. 1 is a graph showing spectral accuracy of a multi- 
domain Chebyshev collocation method, in accordance with one 
5 embodiment of the present invention. 

FIG. 2 is a graph showing a comparison of h-type, p-type, 
and hp-type refinement strategies, in accordance with one 
embodiment of the present invention. 
Q FIG. 3 is a graph showing waveforms obtained using 

HP collocation method with hp-type or p-type refinement, in 
W accordance with one embodiment of the present invention. 

m DESCRIPTION OF THE P RE FERRE D EMBO DIMENTS 

, m — — — —— — — — — — — — i — i I I I ■ ■ ■ 

U 

O An invention is disclosed for a method and device for 

multi-interval Chebyshev collocation for efficient, highly 

J J 5 accurate radio frequency circuit simulation. Numerous 
specific details are set forth in order to provide a thorough 
understanding of the present invention. It will be 

understood, however, to one skilled in the art, that the 
present invention may be practiced without some or all of 
20 these specific details. 
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General Overview Of Circuit Equation Discretization 



Circuit behavior is usually described by a set of N 
nonlinear differential-algebraic equations (DAEs) that can be 
written as 

^» + /(v(0) + «(/) = 0, (1) 

dt 

where Q(v(t)) represents contribution of the dynamic elements 
such as inductors and capacitors, I(v(t)) is the vector 
representing static elements such as resistors, u(t) is the 
vector of inputs, and the vector v(t) contains state variables 
such as node voltages. The periodic steady-state solution is 
the solution of Eq. (1) that also satisfies the two-point 
boundary condition v(T) - v{0) = 0 for some period T. 

We now contrast three ways of discretizing the time- 
derivative operator: linear multistep methods, harmonic 
balance methods, and Chebyshev methods. As representative of 
linear multistep methods we consider the Gear methods. Given 
the values, u{tj) , at the time points: t jf j = 0,1,*",M, to 
perform the discretization at t = t kr we seek an interpolating 
polynomial of degree p based on 

u(tj) , j = k - p, k - p + l,'"r k * l" he time derivative is 
approximated by evaluating the derivative of the interpolating 
polynomial at t = t k . Thus, the time-derivative of the 
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approximate solution is required to match the differential 
equation at the point t k . We call such a point a collocation 
point . 

The Gear methods are implicit in only "one-point-at-a- 
5 time". They are called one-stage methods. At each timepoint, 
for a method of order p, the solution vector at only p + 1 
points is involved, and if p previous timepoints are known, a 
nonlinear system of size N must be solved to find the solution 
at the next timepoint. Each step of a shooting method 
1^ involves the factorization of the circuit Jacobian matrix, 
hi Because this matrix is usually very sparse, it can be factored 
with few fill-ins, in nearly 0(N) time. The order of 
? polynomial approximation p is usually fixed at a tow value, 
IP and increased accuracy is obtained by locally changing the 
ijg timesteps based on estimates of the truncation errors induced 
■ by the order-p polynomial approximation. Note that higher- 
order Gear methods can not be used at points after C° or 
sharp-transition points because approximation of the 
differential operator using a p-order Gear method involves p 
20 preceding points. Higher-order approximation across C° points 
or sharp-transition points leads to oscillation of the 
solution and makes the solution even less accurate than a 
first-order solution. A similar problem occurs if the 
timestep ratios are changed too rapidly. 
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For contrast, consider a particular (time-domain) 
construction of the harmonic balance method, known as Fourier- 
collocation. We seek to form the solution u(t) by an 
interpolating trigonometric polynomial of degree M: 



k=-M 



where I M denotes an M-point interpolation operator. By 
computing the derivative of the trigonometric polynomial 



M-l 

(I M u)'(t)=^iku ,h , (3) 



EL 



at the Fourier collocation points tj, tj = M f j = 0,..., 2M - 1 
iS / and enforcing the differential equation at these points we 
JL have specified the numerical method. We may solve for the 
ftj Fourier coefficients directly, or express the Fourier 

Mb 

h coefficients in terms of the values u(tj) at the Fourier 
collocation points as 

The differentiation procedure may be performed 
efficiently to high order by noting that in the frequency 
domain, the differentiation operator becomes a multiplication 
operator, so by using the FFT to convert the u(tj) to u k , 
20 performing the differentiation, then in-verse transforming, 
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one may accomplish the time-derivative operation in 
0(M log 2 M) operations. Alternatively, one may use this 
prescription to obtain a time-domain differentiation matrix D 
such that DV approximates the value of ^r 1 at the collocation 
5 points, where vector V is a vector of the values of v(t) at 
the collocation points. Then the derivatives can be obtained 
directly by a matrix-vector multiplication DV, requiring O(lrf) 
operations. Note that, in practice, using the FFT does not 
provide much advantage over using the D matrix explicitly 
1(1 unless the order of approximation M is fairly large. 
?® Accuracy in the harmonic balance method is usually 

jl achieved by increasing the order of approximation. For a 
T given set of collocation points, the maximal possible order of 
P approximation is used. With Cf° waveforms, this has the very 
l|f desirable effect of achieving spectral convergence, meaning 
^ that as the spacing between the collocation points decreases, 
asymptotically the error decreases faster than any polynomial 
approximation (in particular, faster than the Gear methods) . 
However, note that although the expansion in Fourier series 
20 conveniently means automatically enforcing the periodic 
boundary condition, it also means that the basis functions, 
sines and cosines, must cover the entire domain. Unlike the 
Gear methods, which couple only local points, every point in 
the solution interval [0, T] is coupled to every other through 
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each basis function. The spectral convergence of the Fourier 
series is very desirable for smooth waveforms and nearly 
linear circuits, but unfortunately, a single discontinuous 
point can destroy the spectral convergence properties. 
Because of the global nature of the basis functions, spectral 
convergence is very difficult to recover once interrupted. 
Most semiconductor device models possess some degree of 
discontinuity in the derivatives, and so true spectral 
convergence is rarely observed in practice. 



Chebyshev Polynomials 

A preferred class of basis functions that also possess 
good approximation properties, including the possibility of 
spectral convergence property, are the Chebyshev polynomials. 
However, Chebyshev polynomials have the advantage that they 
can be applied when the simulation domain has been decomposed 
into multiple intervals in order to resolve complex or 
singular behavior, as has been done for electromagnetic wave 
problems in complex geometries. For more information on 
20 electromagnetic wave problems in complex geometries, see B. 
Yang, D. Gottlieb and J.S. Hesthaven, "Spectral Simulations of 
Electromagnetic Wave Scattering," J. Comput. Phys . , vol. 135, 
pp. 216-230, 1997. 
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The following discussion provides the basics of the 
Chebyshev collocation method. Given the values, u{tj) , at the 
Chebyshev Gauss-Lobatto collocation points, 

tj = cos j = 0,..., M f an interpolating polynomial of degree M 
5 is sought: 



I M u(t)^uJ k (t) 9 (5) 

k=Q 

and compute the derivative of the polynomial as 

M 

{I M u)'(t) = ^u' k T k (t). (6) 



jfc=0 



m The coefficients, u' k , are computed from the coefficients u k by 

lift using the backward recurrence relation. 

9 "m +1 ="' = 0 (7) 

n n 

K c t i/;-i5; +2 + 2(i + l)w M ^U,...,M-l 

For more detail, see B. Fornberg, W A Practical Guide to 
Pseudospectral Methods/' Cambridge University Press, 1996. 

We can also obtain a derivative matrix D by noting that 
15 the interpolating polynomial can be expressed in terms of a 
Lagrange inter-polation polynomial gj(t). The entries of the 
Chebyshev differentiation matrix D are computed by taking the 
analytical derivative of gj{t) and evaluating it at the 
collocation points t* for j, k = 0,..., M, i.e., D k j = g'j{t k ) . 
20 The entries of the matrix are 
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A 



kk 




k*0,M 



(8) 



00 



MM ~~ 



2M 2 + 1 
6 



where c k = 2 when k = 0, M and c* = 1 otherwise. It is pointed 
out in W. S. Don and A. Solomonoff, "Accuracy and Speed in 
Computing the Chebyshev Collocation Derivative, " IMA technical 
report, that calculating the matrix directly in this form 
results in unnecessarily large numerical error, and a 
formulation that yields more accuracy is given. 

An important aspect of the multi-interval Chebyshev (MIC) 
approach is to break the simulation domain [0, T] into 
multiple intervals, [0, Ti] , [T X/ T 2 ] r [ T mr T] . In each 
interval [T if T x+ i] the Chebyshev collocation method is used to 
discretize the differential equation. Continuity of the 
solution is enforced at the boundaries of any two intervals. 
For the first and last interval, this has the effect of 
enforcing the periodic boundary condition. The reason this 
decomposition is useful is that waveforms in circuit 
simulation typically have regions of different properties. 
There are smooth regions where high-order methods can be used 
and there are sharp-transition regions where interval 
boundaries must be placed. Our particular combination of 
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high-order Chebyshev discretizations in multiple time 
intervals provides several important advantages. 

First, as in harmonic balance, high accuracy is 
achievable because in each interval all the desirable 
approximation properties of Chebyshev polynomials can be 
exploited. Chebyshev polynomials are usually close to the 
optimal polynomial approximation of a given order, and 
spectral convergence can be achieved if order-based (p-type) 
refinement is used in each interval. Because of the rapid 
convergence, once a proper set of intervals is chosen high 
accuracy can be obtained at little marginal cost. Conversely, 
there will be few intervals in regions where the solution is 
smooth, because the intervals there can be quite large. 

Second, as in the lower-order finite difference schemes, 
the timestep intervals can be adaptively chosen to resolve 
nonlinear or abrupt transition behavior. The impact of 
nonlinear behavior can some-times be difficult to diagnose. 
Solution waveforms, such as voltages, may be quite smooth, but 
residuals, such as currents, will be much more difficult to 
represent with a Fourier series if the circuit equations are 
very nonlinear. In the MIC method, if the intervals are 
chosen such that discontinuities occur only at the interval 
boundaries, then the discontinuities can be handled 
efficiently. Around a transition or strongly nonlinear 
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region, more intervals can be placed in a denser timestep 
spacing, and possibly the order of the method lowered. 

Third, the method has excellent numerical stability 
properties. It can be shown that the scheme of the present 
5 invention is a member of a particular class of implicit Runge- 
Kutta (IRK) methods, the collocation-IRK methods. For more 
information on solving differential equations, see E. Hairer 
and G. Wanner, "Solving Ordinary Differential Equations II," 
Q Springer-Verlag, 1991. Higher-order implicit Runge-Kutta 
t6 methods are multi-stage methods, that is, they are implicit in 
^ more than one point at a time. Because of the multiple 
^ implicit states, and unlike multi-step methods such as the 
f*i Gear methods, IRK methods can be constructed that are A-stable 
p even at high order. A-stable means a scheme that is 
|Q5 numerically stable for systems with poles located anywhere in 
the left-half-plane, including the entire imaginary axis. The 
particular schemes used in this discussion are not A-stable, 
but it has been proved that the schemes of all orders used 
here are what is called A(a) -stable, meaning that if X is a 
20 pole of the system being integrated, then the numerical scheme 
is stable with timestep h for all z = hX in the left-half- 
plane, except possibly within an angle a of the imaginary 
axis. For the methods of the present invention, a is usually 
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less than one degree, and the stability region only excludes a 
small semi-circle near the imaginary axis, which can be 
avoided by appropriate choice of timesteps. 

In addition, because the Gauss-Lobatto-Chebyshev points 
are preferred for collocation (implying a collocation point is 
always placed at the forward interval boundary) , the 
particular construction here gives methods that are L- or 
stiffly stable, just as the Gear methods are. Accordingly, 
the method is suitable for solving differential-algebraic 
equations as initial value problems. This property allows the 
use of shooting techniques to solve boundary value problems or 
construct preconditioners . In fact, in a certain sense that 
is important in practical applications, the MIC approach is 
more stable than the Gear methods. Since unlike the Gear 
methods, the MIC scheme retains excellent stability properties 
even when the timesteps are rapidly varied, higher-order 
schemes can be used near sharp transition points. Finally, 
the discretization scheme is robust. Unlike harmonic balance 
schemes, where sufficiently sharp transitions can induce 
oscillations in the global basis that are next to impossible 
to eliminate, in the worst possible case the discretization in 
a single MIC interval degenerates to the backward-Euler scheme 
which is robust and well-tested. In practice, the need to 
drop below the second order scheme is rare. 
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Efficient Implementation of Chebyshev Collocation Methods 

For simplicity of presentation, the following discussion 
involves Chebyshev-collocation methods for the case of two 
intervals: I x = [0 , Zi] and I 2 = [Ti,T]. The number of 
collocation points in Ii is L and the number of collocation 
points in I2 is M. The solution in Ii is represented by 

V = Vof Vi, V 1/ - t , V L and the solution in I 2 is represented by 

V = V 0/ Vi,"-, V m . The continuity condition V h = % is enforced 
at the internal interval boundary, and the periodicity 

condition Vb = relates the two interval endpoints. Assume 
that D = (dij)|i,j - 0,l f "',L is the differentiation matrix in 

Ii and B = {dij)\i,j = 0,1, "',M is the differentiation matrix 
in I2, Eq. (1) is discretized into the following equation: 
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d n l - d u l 0 0 d l0 l 



d Ll l *" 
0 0 



d l0 I d n l 



0 <y 



Q + I + U = 0, 



(9) 



/ + u = 0, 

where I is the identity matrix. Newton ! s method is used to 
solve this nonlinear equation. By defining the matrices 



dV dV 



(10) 
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at each collocation point we can write the Jacobian matrix J 



for the Newton iteration as 



0 d w C 0 ^n c i +G i 



0- 0 d M0 C 0 d Ml C i - d MM C M +G M- 



(ID 



Note that C 0 = C M and C L - C Q . 

At each Newton iteration step, we must solve a linear 



system 

JAU = -E 9 



(12) 



where E is the residual of the previous solution. This step 
is problematic because in each interval the discretization of 
the differentiation operator is global. That is, the 

approximation of ^ at every collocation point involves 
q(v(ti) at all the collocation points t±. Consequently, the 
differentiation matrix, like the Fourier collocation 
differentiation matrix, is a full matrix (on each interval) . 
The resulting larger, less-sparse equations systems lead to 
difficulties in large-scale circuit simulation and is the 
reason IRK methods are not widely used to date. Ironically, 
it is the sheer scale of RF simulation problems that makes IRK 
methods more suitable in the RF context. Matrix-implicit 
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Krylov subspace methods , such as Generalized Minimum Residual 
(GMRES) or Quasi-Minimal Resudual (QMR) , are now commonly used 
to simultaneously solve the linear system of equations (Eq. 
11) resulting from discretization of the circuit equations 
over a fundamental period T. The present invention shows that 
it is possible to easily construct a preconditioner that is 
more effective than those typically used in harmonic balance 
and comparable to the preconditioners used for shooting 
methods, thereby assuring rapid convergence of the iterative 
solvers. For example, on a similar test (the DC-DC 

converter) , with a comparable approach to preconditioner 
construction, the method of the present invention requires 
two orders of magnitude fewer GMRES iterations than reported 
by 0. J. Nastov and J. K. White, "Time-Mapped Harmonic 
Balance," Proc. 36th Design Automation Conference, New 
Orleans, LA, June, 1999. 

In each interval a "local" preconditioning approach is 
preferred that is motivated by techniques used for large-scale 
harmonic balance analysis. But because [0 f T] is decomposed 
into multiple intervals, a preconditioner can be devised that 
is much more efficient than the corresponding harmonic balance 
preconditioner. In this way, the "global" preconditioning 
approach is similar to that used in R. Telichevesky, K. 
Kundert and J. K. White, "Efficient AC and Noise Analysis of 
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Two-Tone RF Circuits," Proc. 33rd Design Automation 
Conference, June, 1996, for shooting methods, where a 
preconditioner is obtained for finite difference equations by 
dropping the top-right corner entries in the left-hand-side 
matrix of the linear system. 

In harmonic balance, it is common to construct a 
preconditioner by assuming the capacitance and conductance 
matrices C, G respectively can be well represented by 
(piecewise) constant approximations. Both these techniques 

are adopted in the present invention. For example, a single C 
is used as an approximation of all the capacitance matrices in 

interval Ii and C as an approximation of all the capacitance 
matrices in interval I 2 , the conductance matrices being 
treated similarly. When this is done, note that the two large 
blocks of the preconditioning matrix corresponding to the two 
intervals are only coupled by the subdiagonal block. Hence, 
it suffices to discuss how to invert the top diagonal block. 
To accomplish this inversion method, consider forming the 
Schur decomposition of the relevant portion of the Chebyshev 
differentiation matrix such that 

[dX^^ mU ^ d3) 

where T is a Schur matrix and U is a unitary matrix. T has 2- 
by-2 blocks on the diagonal and all other subdiagonal entries 
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are zeros. By applying this decomposition to the top block of 
J, one can reduce it to an almost upper block-triangular 
structure with 2N-by-2N blocks on the diagonal. These 2N-by- 
2N blocks on the diagonal can be further transformed into an 
N-by-N complex matrix that has the same sparsity structure as 
a single-timepoint circuit Jacobian as generated in a 
transient analysis. It can be factored efficiently, with few 
fill-ins, just as the Jacobian matrix in a low-order method. 

The advantage of this preconditioner over harmonic 
balance preconditioners is that one can approximate the 
capacitance and conductance matrices with averages separately 
in each interval at little computational cost, in contrast to 
using global averages which make the preconditioner very 
inefficient in simulations of nonlinear circuits. Because the 
intervals are chosen to resolve the nonlinear behavior, we can 
expect the cost of applying the preconditioner to be 
relatively independent of the degree of nonlinearity of the 
circuit. If the circuit is very nonlinear, there will be very 
many intervals, yet precisely for this reason the 
preconditioner will still be a good approximation to the 
original matrix, and in fact will be nearly as cheap to apply 
as the shooting-method preconditioner. 

In a practical implementation directly solving the MIC 
equations is not preferred. Instead, a low-order shooting- 
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Newton method is used to obtain an initial approximation to 
the solution waveform. For steady-state problems, the 
shooting-Newton method has a good global convergence rate in 
periodic steady state (PSS) simulations and can rapidly reach 
5 the neighborhood of the steady state. 

From the result of the low-order solution method, the 
time interval [0, T] is decomposed into a number of 
subintervals. The objective of this first decomposition is to 
^ avoid sharp transitions and C° points in each interval. Next, 
l8 the order of accuracy desired in each interval is determined. 
UJ At this point, if an interval contains many time points of the 
5 low-order shooting method, the solution in this interval must 
*' be relatively smooth, and a relatively high order in this 
0^ interval will be chosen. If an interval contains only a few 
l!| time points, it means the interval is in a sharp-transition 
^ region and low-order IRK methods in the interval can be 
chosen. Once the circuit equations have been discretized in 
all the intervals, and the boundaries are connected with 
continuity conditions and the periodicity conditions, then the 
20 Newton iteration can proceed. 

Once an initial solution is obtained, an adaptive 
refinement scheme is used to split the intervals in any 
additional sharp-transition regions that may appear and 
increase the order of approximation in other regions where 
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accuracy is not satisfactory. In this stage, both h-type 
(interval decomposition) and p-type (order increase) 
refinement are used to achieve the required accuracy in the 
most efficient way. 
5 Error control can be performed by examining the trailing 

Chebyshev coefficients as indicators of potential numerical 
truncation error. Initial experiments have been conducted 
using the following simple error criteria, 

O (^ + C N )h s /T<e, (14) 

where N s is the order of the Chebyshev method used in inerval 
I! s, c Na _! and c Ns , are the N s - 1th and N s th th Chebyshev expansion 
coefficients, normalized by the local maximum in all 

G 

P intervals, h s is the length of the interval, T is the time 
5f period, and a is the desired level of accuracy. This takes 
# into account both h-type refinement, h s /T, and the p-type 

C N - 1 • 
refinement, — ^ + C^ s/ which is proven to be effective in 

numerical experiments. If the error criteria is not satisfied 
in a given interval, then we check the absolute value of Cu s -1 
and C Ns . If ICn s -iI > I Cn s I , we increase the order of the 
20 method in the interval, as this indicates that the order of 
Chebyshev coefficients is decreasing, so it is reasonable to 
seek spectral convergence. Otherwise, the solution is likely 
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non-smooth, and we decompose the interval into two new 
intervals at the middle point. 

Numerical Results 

In this section, the spectral accuracy of the multidomain 
Chebyshev collocation method in periodic steady state computa- 
tions is provided. Efficiency of using h-type (domain 
decomposition) and p-type (polynomial order increase) 
refinement in resolving waveforms with sharp transitions is 
shown. The method's ability to resolve sharp transitions with 
high accuracy is verified. The efficiency of the 

preconditioning technique is also demonstrated by the solution 
of the linear system using the matrix-free Krylov-subspace 
approach. 

In the first example, the spectral accuracy of the 
multidomain Chebyshev collocation methods is shown by 
computing the periodic steady-state solution of a crystal 
filter. 

FIG . 1 shows how the method of the present invention can 
achieve machine precision using only 32 time points (we use 2 
intervals and 16th order in each interval), in accordance with 
one embodiment of the present invention. Spectral convergence 
is evident in the exponential decay of the error. 
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In the second example, the efficiency of using both h- 
type and p-type refinement in resolving solutions with sharp 
transitions is provided. Nonlinearity is introduced into the 
previous filter example by placing a mixer before it. Three 
refinement strategies are compared: h-type refinement only, 
p-type refinement only, and hp-type (h-type and p-type) 
refinement . 

FIG. 2 shows how both p-type refinement and hp-type 
refinement can achieve spectral accuracy, in accordance with 
one embodiment of the present invention. In the hp-type 
refinement, Chebyshev polynomials are used with an order 
varying from 2 to 16 in the intervals. The same level of 
accuracy is achieved that is achieved using p-type refinement 
where the order as high as 66 is used in two intervals. Note 
that for the same total number of time steps, the hp-type 
refinement is more efficient. One can see from FIG. 2 that 
both hp-type refinement and p-type refinement is superior to 
h-type refinement of a second order method. Advantages of 
using high-order methods are evident from this example. 

In the last example, the periodic steady-state of a very 
nonlinear DC-DC converter is calculated, similar to a circuit 
analyzed in 0. J. Nastov and J. K. White, M Time-Mapped 
Harmonic Balance," Proc. 36th Design Automation Conference, 
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New Orleans, LA, June, 1999. The final mesh of the proposed 
approach is obtained through three steps of hp-refinement. 

By using both h-type and p-type refinement, the method 
selects high-order approximations, with large intervals, in 
5 the smooth region and low-order approximations, with small 
intervals, in the sharp-transition region. 

FIG. 3 shows computed waveforms obtained using 
collocation method with hp-type or p-type refinement, in 
pi accordance with one embodiment of the present invention. The 

1|| result is oscillation free and very accurate, as noted from 

%j 

Ly the comparisons of the error in the output voltage presented 

, 

® in Table 1 below. 





Method 


# Time Steps 


Relative Error 




2 nd -order Gear 


2500 


5.0372e - 5 




p-type High- 
order 


128 


8.7343e - 2 




hp-type High- 
order (2-16) 


198 


4.8453e - 6 



15 Table 1: DC-DC converter example. Error of 
the DC output voltage for hp-type high- 
order, p-type high-order, and low-order 
method. 



20 The second order Gear method requires two orders of 

magnitude more timepoints to reach comparable accuracy, while 
using global high-order basis functions, such as used in 
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harmonic balance, leads to large errors. These large errors 
can be explained by examining the strictly high-order solution 
shown in FIG. 3, which shows about the best result obtainable 
using only p-type refinement in two intervals. Oscillation 
phenomena are visibly noticeable. 

It is difficult to make precise comparisons of the 
efficiency of the method, because the shooting solution is 
used to achieve the initial mesh refinement, and because the 
cost varies depending on how many steps of mesh refinement are 
needed, the choice of approximation order in each subinterval, 
etc. However, the method of the present invention generally 
only uses 100 - 150% of the computation time of the initial 
low-order shooting-Newton guess to attain an accuracy of 
orders of magnitude higher. The GMRES iterative solver 
requires somewhat more GMRES iterations (usually between about 
1 and 3 factors) per Newton iteration than the low-order 
shooting-Newton method. However, few Newton iterations beyond 
the shooting stage are typically required to reach 
convergence . 

In the DC-DC converter example, the number of GMRES 
iterations required is only around 10, compared with hundreds 
or even thousands of iterations required by the harmonic 
balance method in 0. J. Nastov and J. K. White, ^Time-Mapped 
Harmonic Balance," Proc. 36th Design Automation Conference, 
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New Orleans, LA, June, 1999. More importantly, because many 
fewer timepoints are used to achieve high accuracy than other 
methods, the memory usage of our scheme, which is usually the 
limiting factor in determining capacity of an RF simulator, is 
potentially very much lower. As a study of efficiency on a 
larger example, what is used is the benchmark receiver circuit 
from R. Telichevesky, K. Kundert and J. K. White, "Efficient 
AC and Noise Analysis of Two-Tone RF Circuits/' Proc. 33rd 
Design Automation Conference, June, 1996. 

With the accuracy requirement at a fixed high precision 
(four digits in the fourth harmonic), the MIC method has shown 
a performance improvement over the low-order shooting-based 
method of 13X in CPU time and 20X in timepoint usage 
(proportional to memory) . If instead the number of timepoints 
is fixed at a similar number in both methods, an additional 
150% of the time of an initial low-accuracy shooting solution 
may be needed by the MIC method to increase the accuracy by 
300X. 
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System And Method Implementation 

Portions of the present invention may be conveniently 
implemented using a conventional general purpose or a 
specialized digital computer or microprocessor programmed 
according to the teachings of the present disclosure, as will be 
apparent to those skilled in the computer art. 

Appropriate software coding can readily be prepared by 
skilled programmers based on the teachings of the present 
disclosure, as will be apparent to those skilled in the software 
art. The invention may also be implemented by the preparation 
of application specific integrated circuits (ASIC's) or by 
interconnecting an appropriate network of conventional component 
circuits, as will be readily apparent to those skilled in the 
art . 

The present invention includes a computer program product 
which is a storage medium (media) having instructions stored 
thereon/in which can be used to control, or cause, a computer to 
perform any of the processes of the present invention. The 
storage medium can include, but is not limited to, any type of 
disk including floppy disks, mini disks (MD's), optical discs, 
DVD, CD-ROMS, micro-drive, and magneto-optical disks, ROMs, 
RAMs, EPROMs, EEPROMs, DRAMs, VRAMs, flash memory devices 
(including flash cards) , magnetic or optical cards, nanosystems 
(including molecular memory ICs), RAID devices, remote data 
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storage/archive/warehousing, or any type of media or device 
suitable for storing instructions and/or data. 

Stored on any one of the computer readable medium (media) , 
the present invention includes software for controlling both the 
hardware of the general purpose/specialized computer or 
microprocessor, and for enabling the computer or microprocessor 
to interact with a human user or other mechanism utilizing the 
results of the present invention. Such software may include, 
but is not limited to, device drivers, operating systems, and 
user applications. Ultimately, such computer readable media 
further includes software for performing the present invention, 
as described above. 

Included in the programming (software) of the 
general/specialized computer or microprocessor are software 
modules for implementing the teachings of the present invention, 
including, but not limited to, applying Chebyshev collocation to 
boundary-value differential equations to discretize a set of BVD 
equations, and forming a solution to the set of BVD equations 
based on the discretized BVD equation, according to processes of 
the present invention. 

In the foregoing specification, the invention has been 
described with reference to specific embodiments thereof. It 
will, however, be evident that various modifications and 
changes may be made thereto without departing from the broader 
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spirit and scope of the invention. The specification and 
drawings are, accordingly, to be regarded in an illustrative 
rather than a restrictive sense. 
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